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1.0  INTRODUCTION 


The  subject  of  digital  filters  Is  somewhat  confused  because  a digital 
filter  Is  usually  considered  to  be  a computational  procedure  rather  than 
a physically  recognizable  piece  of  hardware.  These  computational  pro- 
cedures vary  greatly  in  complexity  and  purpose.  Some  may  digitally  sim- 
ulate a real  electrical  or  mechanical  filter.  Others  may  simulate  ideal- 
ized filters  which  could  not  be  easily  made  into  hardware,  and  still 
others  may  be  simple  numerical  smoothing  techniques  not  related  to  any 
actual  or  ideal  filter. 

The  methods  presented  in  this  paper  are  based  on  the  use  of  the 
reversible  time  to  frequency  transformation  - the  'Fast  Fourier  Transform' 
(FFT) . To  put  this  method  into  proper  perspective,  some  traditional 
methods  are  briefly  discussed  in  the  next  section. 
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2.0  TRADITIONAL  METHODS 


Various  numerical  procedures  for  smoothing  discrete  data  traces  have 
been  called  digital  filters.  Some  of  these  methods  simply  replace  each 
data  point  with  a weighted  average  of  the  surrounding  points.  The  com- 
putations are  usually  done  on  a digital  computer.  Another  method  treats 
filtering  as  an  initial  value  problem.  That  is,  a digital  simulation  of 
a real  mechanical  or  electrical  system  is  used  as  a filter.  The  discrete 
data  record  to  be  filtered  is  considered  to  be  a forcing  function  or 
applied  voltage,  and  a computed  force  or  voltage  at  some  point  in  the 
system  is  the  filtered  data.  A digital  computer  numerically  Integrates 
the  system  differential  equations  to  obtain  the  output. 

These  methods  have  several  disadvantages.  The  actual  effect  upon 
the  frequency  content  of  the  data  being  filtered  is  not  intuitively 
obvious.  The  effective  transfer  function  of  the  simple  averaging  process 
depends  upon  the  sample  spacing  as  well  as  the  particular  set  of  weighting 
factors  used.  In  the  numerical  integration  method,  the  parameters  of  the 
system  must  be  carefully  chosen  to  obtain  the  desired  frequency  response 
and  scale  factors.  Both  methods  may  produce  some  undesirable  amplification 
and/or  shifts  in  the  frequencies  in  the  pass  band,  and  the  desired 
transfer  function  can  usually  only  be  approximated.  Also,  there  may  be 
several  pass  bands  because  of  the  phenomenon  of  frequency  aliasing.  To 
avoid  thase  limitations,  the  following  FFT  method  is  becoming  popular. 
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3.0  THE  FFT  METHOD 


The  Discrete  Fourier  Transform  (DFT)  Is  the  basic  tool  which  facilitates 
the  methodology  presented  in  this  paper.  The  DFT  is  defined  by: 

N-l 

Ar  - £ \ EXP  (-2ir  jrk/N)  (1) 


for  r ■ 0t  ...  N - 1 


The  usefulness  of  the  DFT  was  limited  because  the  computational  time  was 
too  great , even  on  the  faster  digital  computers.  A breakthrough  occurred 
in  the  mid  60's  with  the  development  of  a much  faster  algorithm  for  com- 
puting the  DFT.  This  method  is  commonly  called  the  'Fast  Fourier  Transform' 
(FFT).*  The  FFT  provides  an  efficient  procedure  for  transforming  a discrete 
time  series  into  a complex  function  of  frequency.  No  information  is  lost 
in  this  transformation;  in  fact,  the  transformation  is  reversible.  This 
reversibility  feature  of  the  FFT  suggests  the  possibility  of  modifying  a 
time  series  by  operating  on  its  FFT  followed  by  an  inverse  FFT.  Filtering 
may  be  accomplished  by  the  convolution  of  the  FFT  and  the  frequency  transfer 
function  of  a filter  followed  by  an  inverse  FFT.  This  is  equivalent  to 
solving  the  convolution  integral  to  obtain  the  filter  output: 


C(t) 


X(t)  h(t-t)  dt 


(2) 


where: 

C (t)  is  the  output 
X (t)  is  the  input 

h (t)  is  the  response  of  the  filter  to  a step  input 

The  direct  solution  of  equation  (1)  is  completely  impractical  if  it 
must  be  done  by  numerical  Integration.  Even  if  the  input  function  x(t) 
is  zero  for  t < o so  that  the  lower  limit  on  the  integral  becomes  zero, 
the  computation  time  is  much  too  great. 


* The  FFT  in  this  study  was  computed  using  subroutine  F0UR1,  written  by 
Norman  Brenner  of  the  MIT  Lincon  Laboratory,  July  1967.  A listing  of 
FOUR!  is  included  in  the  appendix. 
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4.0  IDEALIZED  FILTERS 


Ideal  or  perfect  filtering  may  be  accomplished  by  setting  the  FFT 
coefficients  corresponding  to  frequencies  outside  the  pass  band  to  zero 
followed  by  an  inverse  FFT.  The  passed  frequencies  are  essentially 
unchanged  in  both  magnitude  and  phase.  Subroutines  for  this  type  of 
filtering  (PLPASS,  PHPASS,  and  PBPASS)  which  are  Perfect  Low  PASS,  Perfect 
High  PASS,  and  Perfect  Band  PASS  filters  respectively  have  been  written 
and  are  included  in  the  appendix. 

Examples  of  using  these  filters  follow.  Figure  1 shows  the  test  data 
used  to  demonstrate  each  of  the  filters  presented  in  this  paper.  Figure  2 
shows  the  power  spectrum  of  the  test  data  which  consisted  of  discrete  fre- 
quencies at  1.0,  2.25,  5.0,  10.0,  20.0,  and  30.0  HZ.  This  hypothetical 
data  was  used  because  it  does  contain  a wide  range  of  strong  frequencies, 
thus,  providing  the  opportunity  of  demonstrating  low-pass,  high-pass,  and 
band-pass  filters  all  using  the  same  input  data.  Figure  3 shows  the  out- 
put from  the  Perfect  Low  PASS  (PLPASS)  filter.  The  cutoff  of  3.625  HZ  only 
allows  the  frequencies  at  1.0  and  2.25  to  pass.  Figure  4 shows  the  output 
from  the  Perfect  High-Pass  (PHPASS)  filter  with  a cutoff  frequency  at  15.0 
HZ  allowing  only  the  frequencies  at  20.0  and  30.0  HZ  through.  Figure  5 
shows  the  use  of  the  Perfect  Band  PASS  (PBPASS)  filter.  Note  that  the  low 
cutoff  of  7.5  and  the  high  cutoff  of  25.0  HZ  allows  only  the  frequencies 
at  10.0  and  20,0  HZ  through. 

The  selection  of  cutoff  points  for  filtering  real  data  should  be 
based  on  some  knowledge  of  the  system  from  which  the  data  was  obtained. 

For  example  the  hull  roll  motion  of  a tank  has  a natural  frequency  of 
about  1.0  HZ,  thus,  any  frequency  content  greater  than  5.0  to  10.0  HZ  is 
probably  noise  and  may  be  filtered  out.  One  source  of  noise  in  tank  data 
which  has  been  observed,  is  caused  by  the  track  shoes  contacting  the  ground. 
This  shows  up  on  a power  spectrum  plot  as  a sharp  peak  at  15.0  to  20.0  HZ. 
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Figure  1 Plot  of  Original  Test  Data 


-> 

I 

u> 


o 

o 

• 

o 

ID_ 

CSI 


o 

o 

o 

o. 

CM 


SPECTRAL  PLOT  OF  TEST  DATA 

(UNSMOOTHED) 


UJo 


cuts 


a: 


o2* 

llj 

CL 

CQ 

o 

o 

• 

o. 

LO 


o 

o 


*^.00 


4 


5.00  10.00  lb .00  20.00 

FREQUENCY  (HERTZ) 


25.00 


30.00 


lb. 


00 


Figure  2 Spectral  Plot  of  Test  Data  (Unsmoothed) 


PLOT  OF  DflTfi  FILTERED 


Figure  3 Plot  of  Data  Filtered  with  PLPASS,  BF-3.625 


a 

to 


Figure  4 Plot  of  Data  Filtered  with  PHPASS,  BF-15.0 
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Figure  5 Plot  of  Data  Filtered  with  PBPASS,  FL-7.5,  FH-25.0 


5.0  DIGITAL  SIMULATION  OF  REAL  FILTERS 


Although  the  perfect  filters  are  very  effective  in  removing  unwanted 
frequency  content  of  a signal,  it  may  be  desirable  to  only  attenuate  cer- 
tain frequency  bands  rather  than  suppress  them  completely.  This  will 
retain  more  of  the  features  of  the  original  signal.  This  is  precisely 
what  a real  RC  filter  does.  The  following  is  a simple  RC  filter: 

R 


Neglecting  loading  effects,  the  transfer  function  can  be  easily  obtained: 

F “ 1 + jw/BF 

where  u ■ frequency 

3 - V-i 


BF  - 1/RC 


Plots  of  this  transfer  function  for  two  break  frequencies  are  included  in 
the  appendix.  To  simulate  this  filter,  the  FFT  of  the  input  is  multiplied, 
point  by  point,  by  the  transfer  function  of  the  filter  and  the  inverse  FFT 
is  obtained.  A subroutine  for  this  filter  (RCN1)  has  been  written  and  is 
included  in  the  appendix.  Figure  6 shows  the  output  of  this  filter  for  a 
particular  break  frequency  (BF)  using  the  same  input  data  as  before.  Note 
that  the  high  frequencies  are  still  present,  although  attenuated.  This 
fact  suggests  the  possibility  of  reversing  the  process.  This  may  be 
accomplished  by  the  same  procedure  except  with  the  inverse  transfer  function: 


F"1  (Jw)  - 7^-  - 1 + jw/BF  (4) 

The  computations  are  basically  the  same,  thus,  subroutine  RCN1  has  been 
written  so  that  it  can  also  reverse  the  process.  This  has  been  done  using 
the  data  in  Figure  6;  the  result  shown  in  Figure  7 is  not  distinguishable 
from  the  original  plot.  This  feature  may  be  used  to  recover  data  which 
was  filtered  at  too  low  a frequency  during  recording  or  digitizing. 
Although  the  original  data  could  be  recovered  by  electrical  processing  of 
the  original  analog  data  and  redigitizing,  the  analog  tape  and/or  the 
special  equipment  and  expertise  to  do  this  may  not  be  readily  available. 
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PLOT  OF  DATA  FILTERED 


Figure  6 Plot  of  Data  Filtered  with  RCN1,  BF-5.0 


PLOT  OF  DATA 


Figure  7 Plot  of  Data  Reconstructed  with  RCN1 


Thus,  this  digital  approach  will  usually  be  more  expeditious.  One  potential 
problem  in  trying  to  reconstruct  the  original  signal  in  this  manner  is  that 
any  noise  in  the  signal  is  also  amplified.  To  help  avoid  this  problem  RCN1 
has  an  optional  sharp  cut  off  frequency,  CF.  Any  frequency  components 
greater  than  CF  are  set  equal  to  zero  as  was  done  in  the  perfect  filters 
discussed  before.  If  CF  is  equal  to  or  greater  than  the  Nyqulst  frequency 

( */ 2 * DT) , which  is  the  highest  frequency  that  can  be  considered,  none 
of  the  frequencies  are  cut  off. 

Additional  subroutines  have  been  written  (RCN2  and  RCN3)  and  are 
included  in  the  appendix.  These  subroutines  correspond  to  the  following 
RC  networks: 

C R1 


I XI 

RCN2  RCN3 

RCN2  has  the  following  transfer  function: 


F<J“>  - ir^r: 

where : 


(5) 


u - frequency 


BF  - 1/RC 

J -V1!"" 


This  is  basically  a high-pass  filter;  it,  obviously,  will  not  pass  DC. 
Plots  of  the  transfer  function  for  two  values  of  BF  are  included  in  the 
appendix.  Figure  8 shows  the  result  of  using  RCN2  on  the  test  data 
plotted  in  Figure  1. 

RCN3  has  the  following  transfer  function: 


F(J«) 


1_  BF  + Iid 
A BF/A  + jio 


where: 


A - 1 + R1/R2 
BF  - 1/ (R2  * C) 

j 

id  » Frequency 


(6) 
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Figure  8 Plot  of  Data  Filtered  with  RCN2,  BF-10.0 


r 


This  filter  i9  a low-pass  filter,  similar  to  RCN1;  however,  the  attenuation 
curve  in  the  low  frequency  range  is  much  steeper  and  the  attenuation  of 
high  frequencies  approaches  a fixed  value  of  l/A.  A plot  of  the  transfer 
function,  for  particular  values  of  BF  and  A,  is  included  in  the  appendix. 
The  use  of  RCN3  on  the  test  data  of  Figure  1 is  shown  in  Figure  9.  The 
outputs  of  both  RCN2  and  RCN3  were  successfully  used  to  reconstruct  the 
original  data  by  using  RCN2  and  RCN3  with  their  inverse  transfer  functions. 
The  resulting  plots  were  indistinguishable  from  Figure  1 and  are  shown  in 
Figures  10  and  11  respectively. 
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PLOT  OF  DfiTfl  FILTERED 


Figure  9 Plot  of  Data  Filtered  with  RCN3,  BF-5.0,  A- 10.0 


PLOT  OF  DATA 


Figure  10  Plot  of  Data  Reconstructed  with  RCN2 


PLOT  OF  DATA 

RECONSTRUCTED  WITH  RCN3 


Figure  11  Plot  of  Data  Reconstructed  with  RCN3 


6.0  OTHER  NETWORKS 


Subroutines  for  the  digital  simulation  of  other  RC  networks  could  be 
written  by  a relatively  simple  modification  of  RCN1,  RCN2,  or  RCN3.  The 
procedure  is  to  obtain  the  frequency  transfer  function;  rationalize  the 
denominator;  and,  define  the  transfer  function  by  its  real  and  imaginary 
parts  (AMPR  and  AMPI).  Then  do  the  same  for  the  inverse  transfer  function 
and  simply  replace  the  4 statements  where  AMPR  and  AMPI  are  defined  in  one 
of  the  subroutines  RCN1,  RCN2,  or  RCN3.  Note  that  complex  numbers  in  these 
subroutines  are  treated  as  two  real  numbers  so  they  can  be  used  on  computer 
systems  which  do  not  have  a complex  number  capability.  Note,  also,  that 
in  each  of  the  subroutines,  use  is  made  of  the  fact  that  only  the  imaginary 
part  of  the  transfer  function  changes  sign  for  negative  frequencies.  Thus, 
the  Fourier  coefficients  corresponding  to  negative  frequencies  in  the 
modified  FFT  are  the  complex  conjugates  of  the  coefficients  corresponding 
to  positive  frequencies. 
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7.0  PRACTICAL  APPLICATIONS 


Mathematical  models  of  a wide  range  of  weapon  systems  have  been 
developed.  These  systems  range  from  small  arms  to  complete  tanks  and 
helicopter  gun  ships.  These  models  are  useful  in  sensitivity  studies  and 
in  the  evaluation  of  competing  subsystems.  It  is  customary  to  first  validate 
a math  model  by  comparing  the  output  of  the  model  with  experimental  data 
from  laboratory  or  field  tests  and  refining  the  model  until  satisfactory 
agreement  is  obtained. 

i 

As  the  system  becomes  larger  and  more  complicated  the  accuracy  of 
the  modeling  may  become  more  approximate  because  it  is  not  practical 
or  necessary  to  Include  every  detail  of  a large  system  in  a model. 

For  example,  a model  of  a tank  (HITPRO)  does  not  Include  any  provision 
to  account  for  vibration  from  the  engine  and  drive  train,  or  structural 
vibration  of  the  hull  and  turret  which  are  modeled  as  rigid  bodies. 

This  is  not  considered  to  be  a limitation  of  the  model  because,  for  its 
Intended  purpose,  it  is  not  necessary  to  predict  high  frequency  response. 
However,  high  frequency  vibration  from  these  sources  does  show  up  in  the 
experimental  data.  Also,  experimental  data  may  contain  noise  caused  by 
the  instrumentation,  signal  processing,  recording,  etc.  Therefore, 
when  comparing  the  experimental  data  with  the  model  output  it  is  reason- 
able to  filter  out  these  higher  frequencies.  Figure  12  shows  the  out- 
put from  an  accelerometer  mounted  in  the  turret  of  a MICV  (Mechanized 
Infantry  Combat  Vehicle).  The  spectral  plot  of  this  data,  Figure  13, 
shows  a considerable  amount  of  frequency  content  in  the  range  above  20 
Hertz  which  is  obviously  structural  vibration  and  noise.  Figure  14  shows 
the  result  of  filtering  this  data  with  PLPASS  using  a cut  off  frequency 
of  8 Hertz  so  the  data  could  be  compared  with  the  model  output. 

Figure  15  shows  the  vertical  gun  pointing  error  of  an  M60A2  tank  with 
a stabilized  gun  as  predicted  by  a mathematical  model  of  the  system 
(HITPRO).  It  was  desirable  to  obtain  an  RMS  value  for  the  pointing  error 
as  a measure  of  the  performance  of  the  gun  stabilization  system.  This 
value  was  to  be  compared  with  a similar  value  for  a competing  system. 

There  is  reason  to  believe  that  the  upward  drift  shown  in  Figure  15  was 
caused  by  some  inaccuracy  in  the  model  of  the  system  and  the  inability 
of  the  gunner  to  respond  to  small  errors  - not  the  stabilization  system. 
Therefore,  to  obtain  a meaningful  RMS  value  of  pointing  error  as  a measure 
of  the  performance  of  the  stabilization  system,  it  was  desirable  to  remove 
this  drift  before  computing  the  RMS  value.  If  this  drift  was  linear,  one 
could  simply  subtract  an  appropriate  ramp  function  from  the  data.  The 
drift,  however,  appears  to  be  nonlinear  and,  therefore,  some  other  method 
was  required.  Filtering  out  the  very  low  frequencies  was  found  to  be  an 
effective  method  of  removing  drift  from  this  data.  The  result  of  filtering 
the  data  with  the  Perfect  Band  PASS  (PBPASS)  filter  is  shown  in  Figure  16. 

A low  frequency  cut  off  of  0.25  HZ  was  found  to  be  adequate.  A high  fre- 
quency cut  off  of  25.0  HZ  was  used  because  the  actual  system  being  modeled 
is  not  capable  of  responding  to  frequencies  this  high. 
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Figure  12  Turret  Vertical  Acceleration  (Experimental  Data) 


Figure  13  Turret  Vertical  Acceleration  (Unstnoothed  Spectral  Plot) 


♦ i 


Figure  14  Turret  Vertical  Acceleration  Filtered  with  PLPASS, 


CF-8.0 


TRNK  GUN  POINTING  ERROR  (ELEVATION) 


Figure  15  Tank  Gun  Pointing  Error  (Elevation)  (Data  From  Digital  Simulation) 


TANK  GUN  POINTING  ERROR  (ELEVATION) 


Figure  16  Tank  Gun  Pointing  Error  (Elevation)  Filtered  with  PBPASS,  FL-0.25,  FH-25.0 


8.0  CONCLUSION 


The  techniques  presented  In  this  paper  provide  an  efficient  means  of 
altering  the  frequency  content  of  discrete  data  traces  in  a completely 
predictable  manner.  These  techniques  are  applicable  In  many  diverse 
fields  and  data  processing  requirements.  In  addition  to  the  applications 
to  data  from  combat  vehicles,  these  techniques  have  been  successfully 
used  to  smooth  aircraft  flight  path  data,  as  measured  by  tracking  radar. 
Smoothing  was  necessary  because  this  data  was  used  as  an  input  into  a 
mathematical  model,  and  it  was  physically  impossible  for  an  aircraft  to 
follow  the  flight  path  indicated  by  the  raw  data.  These  techniques  have 
also  been  used  in  the  reduction  of  helicopter  flight  test  data. 
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SUBROUTINE  FOUR1  ( DATA .N . ISION) 

C THE  COOLEY-TUKEY  FAST  FOURIER  TRRN3F0RH  IN  U3R3I  BR3IC  FORTRAN. 

C NOTE—  IT  SHOULD  NOT  BE  NECES3RRY  TO  CHRNOE  RNY  STATEMENT  IN  THIS 

C PROORAH  SO  LONO  RS  THE  FORTRRN  COHPILER  USEO  STORES  REAL  RNO 
C IHRDINARY  PARTS  ADJACENTLY  IN  3TORAOE • 

C TRANSFORH(K)  = SUM( OATA( J )»EXP( ISI0N»2»PI»$QRT( -l )■( J-l )»(K-l ) 

C /N)).  SUMMED  OVER  ALL  J RNO  K FROM  1 TO  N.  ORTA  IS  R ONE- 
C DIMENSIONAL  COMPLEX  ARRAY  (I.E..  THE  REAL  AND  IMAOINARY  PARTS  ARE 

C ADJACENT  IN  STORAOE.  SUCH  AS  FORTRAN  IV  PLACES  THEM)  WHOSE  LENGTH 

C N=2«K.  K .DC .0  (IF  NECESSARY.  APPEND  ZEROES  TO  THE  DATA).  I3I0N 

C IS  -1  OR  -l-  IF  A -1  TRANSFORM  IS  FOLLOWED  BY  A ♦ ! ONE  (OR  VICE 

C VERSA)  THE>  ORIGINAL  ORTA  REAPPEAR.  MULTIPLIED  BY  N.  TRANSFORM 

C VALUES  ARE  RETURNED  IN  ARRAY  ORTA.  REPLACING  THE  INPUT.  THE  TIME 

C IS  PROPORTIONAL  TO  N»L002(N).  RATHER  THAN  THE  NAIVE  N«2. 

C ACCURACY  IS  ALSO  GREATLY  IMPROVED.  THE  RMS  RELATIVE  ERROR  BEING 
C BOUNDED  BY  6«3QRT( 2 )»L002( N )«2aa( -B  ) • WHERE  B IS  THE  NUMBER  OF 

C BITS  IN  THE  FLOATING  POINT  FRACTION.  WRITTEN  BY  NORMAN  BRENNER  OF 

C MIT  LINCOLN  LABORATORY.  JULY  1967.  THIS  IS  THE  SHORTEST  VERSION 

C OF  THE  FFT  KNOWN  TO  THE  AUTHOR.  FASTER  PROGRAMS  F0UR2  AND  FOURT 

C EXIST  THAT  OPERATE  ON  ARBITRARILY  SIZED  MULTIDIMENSIONAL  ARRAYS. 

C SEE—  IEEE  AUDIO  TRANSACTIONS  (JUNE  1067).  SPECIAL  ISSUE  ON  FFT. 

DIMENSION  OATA(l) 

IP0=2 

IP3=IP0aN 

!3REV=l 

00  60  13=1. IPS. IPO 
I F ( I3-I3REV 110.20.20 
10  TEMPR=OATA( 13) 

TEMPI=DATA(I3M) 

OATAC 13  )=OATA( I3REV ) 

0ATA(I3*i)=0ATA(I3REV*l) 

0ATA( I3REV )=TEMPR 
0ATA(I3REV*1)=TEHPI 
20  IPt=IP3/2 

30  IF( I3REV-IPI ) 50 .60 .40 
40  I3REV= I3REV-IPI 
IPl=!Pl/2 

I F ( IPl-IPO 150.30 .30 
60  I3REV=I3REV*IP1 
IPl=IPO 

60  IF(IP1-IP3 170.100 .100 
70  IP2=IP1»2 

THETA=6  »2B3 186307/FLOAT ( I SION* IP2/IP0) 

SINTH=SIN( THETA/2.) 

W3TPR=-2.»3INTH»3INTH 
W3TP 1=3 IN( THETA) 

WR=1 . 

WIsQ  • 

00  90  11=1. IPI. IPO 
00  BO  13=11 • IPS. IP2 
I2A=I3 
I2B=I2AWPl 

TEMPR=WR»OATA( I2B)-WI«0ATA( I2B*i ) 

TEMPI=WRaOATA( I2B*l )*MIaOATA( I2B) 

OATA( I2B )=OATA( I2A)-TEMPR 
OATA( I2B>t )=OATA( I2A*l l-TEMPI 
OATA( I2A)=0ATA( I2A)*TEMPR 
BO  0ATA(I2AM)=0ATA(I2A*1)*TEHPI 
TEMPR=WR 

WR=WR»WSTPR-WI»WSTPI ♦WR 
90  WI=WIaWSTPR*TEMPRaWSTPI*MI 
IPUIP2 
00  TO  60 
100  RETURN 
ENO 
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SUBROUTINE  PLPASS(X.N.OT.CF) 


LOW-PASS  OIOITfll  FILTER— U3E3  THE  COOLEY-TUKEY  FfiST  FOURIER  TRANSFORM 
(F0UR1 ) • X 13  THE  ORTA  TRACE  TO  BE  FILTEREO  AND  IS  ASSUMED  TO  OCCUPY  THE 
FIRST  HALF  OF  THE  X ARRAY.  IT  IS  NOT  NECESSARY  TO  INITIALIZE  THE  SECOND 
HALF  OF  THE  X ARRAY  OR  TO  PUT  IT  INTO  COMPLEX  FORM  AS  REQUIRED  BY  FOURl. 

X MUST  BE  DIMENSIONED  2«N  IN  THE  CALLING)  PROORAM  WHERE  N IS  THE  NUMBER  OF 
OATA  POINTS.  NOTE  N MUST  EQUAL  2*«K  WHERE  K IS  AN  INTEOER>ZERO . (APPEND 
ZEROS  TO  THE  OATA  IF  NECESSARY.) 

THE  FILTEREO  OATR  IS  RETURNED  TO  THE  FIRST  HALF  OF  THE  X ARRAY. 

OT  IS  THE  SAMPLINO  RATE  Itf  SECONDS. 

CF  IS  THE  CUTOFF  POINT  (HZ).  OTHER  ROUTINES  EXIST  (PMPASS  ANO  PBPASS) 
WHICH  ARE  HIOH-PASS  ANO  BAND-PASS  FILTERS  RESPECTIVELY. 

OIMEN3 ION  X(l) 

NQc2"N 

FW=FLOAT(N) 

00  1 1=1. N 

X(N0-2«I*l)=X(N-I*l) 

1 X(N0-2*I*2)=0.0 
CALL  FOURKX.N.l) 

IF=CF»FN»OT 

It=IF*2 

I2=N-IF 
00  2 1=11.12 
X( 2« 1-1 )=0 .0 

2 X(2«I)=0.0 

CALL  FOURl(X.N.-l) 

00  3 1=1. N 

3 X(I)=X(2«I-n/FN 
RETURN 

END 
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SUBROUTINE  PHPRSSCX.N.DT .CF) 


HIOH-PRSS  0I0ITRL  FILTER— USES  THE  COOLEY-TUKEY  FRST  FOURIER  TRRNSFORM 
(FOUR!).  X IS  THE  ORTR  TRRCE  TO  BE  FILTEREO  RNO  IS  RSSUMEO  TO  OCCUPT  THE 
FIRST  HRLF  OF  THE  X RRRRY . IT  IS  NOT  NECESSRRY  TO  INITIRLIZE  THE  SECOND 
HRLF  OF  THE  X RRRRY  OR  TO  PUT  IT  INTO  COMPLEX  FORM  RS  REOUIREO  BY  FOUR1. 

X MUST  BE  OIMENSIONEO  2»N  IN  THE  CRLLINO  PROORRM  WHERE  N IS  THE  NUMBER  OF 
ORTR  POINTS.  NOTE  N MUST  EQURL  2»»K  WHERE  K IS  RN  INTEOER>ZERQ . (RPPENO 
ZEROS  TO  THE  ORTR  IF  NECESSRRY.) 

THE  FILTEREO  ORTR  IS  RETURNEO  TO  THE  FIRST  HRLF  OF  THE  X RRRRY. 

OT  IS  THE  SRMPLINO  RRTE  Itf  SECONDS. 

CF  IS  THE  CUTOFF  POINT  (HZ).  OTHER  ROUTINES  EXIST  (PLPRSS  RNO  PBPR3S ) 
WHICH  RRE  L0W-PRS3  RNO  BRND-PRS3  FILTERS  RESPECTIVELY. 

OIMENS ION  X(l) 

N0=2»N 

FN=FLORTCN) 

00  I 1=1. N 
X(N0-2»IM)=X(N-I*1) 

1 XCN0-2»I*2)=0.0 
CALL  FOUR  1 ( X .N . 1 ) 

IF=CF  »FN*OT 
I2=N-IF*-2 

00  2 1 = 1 . IF 
X( 2» 1-1  )=0.0 

2 X( 2» I )=0 .0 
00  3 1=12. N 
X(2»I-l)=0.0 

3 X(2»I)=0.0 

CALL  FOURKX.N.-l) 

00  4 1=1. N 

4 X(1)=XC2»I-1)/FN 
RETURN 

END 
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SUBROUTINE  PBPASS(X.N.OT.FL.FH) 


BRN0-PRS3  OIOITRL  FILTER— USES  THE  COOLEY-TUKEY  FRST  FOURIER  TRANSFORM 
CF0UR1).  X IS  THE  ORTR  TRACE  TO  BE  FILTEREO  RNO  IS  ASSURED  TO  OCCUPT  THE 
FIRST  HALF  OF  THE  X ARRAT.  IT  IS  NOT  NECESSARY  TO  INITIALIZE  THE  SECOND 
HALF  OF  THE  X ARRAY  OR  TO  PUT  IT  INTO  COHPLEX  FORH  AS  REOUIREO  BY  FOUR!. 

X MUST  BE  DIMENSIONED  2«N  IN  THE  CALLING  PROGRAM  NHERE  N IS  THE  NUMBER  OF 
DATA  POINTS.  NOTE  N MUST  EQUAL  2»»K  WHERE  K IS  AN  INTEOER>ZERO.  CAPPENO 
2EROS  TO  THE  DATA  IF  NECES3RRY . ) 

THE  FILTERED  DATA  IS  RETURNED  TO  THE  FIRST  HALF  OF  THE  X ARRAY. 

OT  IS  THE  SAMPLINO  RATE  irf  SECONDS. 

FL  IS  THE  LOW  FREQUENCY  CUTOFF  POINT  (HZ).  AND  FH  IS  THE  HIGH  FREQUENCY 
CUTOFF  POINT.  OTHER  ROUTINES  EXIST  (PLPA33  ANO  PHPAS3)  WHICH  ARE 
LOW-PASS  AND  HIGH-PASS  FILTERS  RESPECTIVELY. 

DIMENSION  X(l) 

N0=2»N 

FNrFLOAT(N) 

00  1 1=1. N 

X(N0-2»I*n=X(N-I*l) 

1 X(N0-2»I*2)=0.0 
CALL  FOUR  1 (X »N . 1 ) 

IF=FL»FN»OT 

I2=N-IF«-2 

00  2 1=1. IF 
X( 2»I-l )=0 .0 

2 X(2»n=0.0 
00  3 1=12. N 

x(2»i-n=o.o 

3 X(2«I)=0.0 
IF=FH»FN»OT 
I 1=IF*2 
I2=N-IF 

00  4 1=11.12 
XC2»I-l)=0.0 

4 X(2»I)=0.0 

CALL  F0UR1CX.N.-1) 

00  G 1=1. N 
6 X(  n=X(2«I-l)/FM 
RETURN 
ENO 
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SUBROUTINE  RCNHX.N.OT.BF .CF. I3I0N) 

C 

C R SIMPLE  L0W-PR33  RC  NETWORK  R3  SHOWN  BELOW t 


C 

C 

C R 

C 0 VVVVV-* 0 FREQUENCY  TRRN3FER 

C \ FUNCTION 

C \ 

C INPUT  OUTPUT  F( JW)=( l/( l^JW/BF) )»»I3ION 

C C 

C \ WHERE! 

C \ BF=l/(R»C) 

C 0 ♦ 0 J=30RTC-1) 

C W=FREQUENCY  (HZ) 

C 

c 

C 13 ION  13  *1  OR  -l i IF  13 ION  13  *1  THE  X VECTOR  13  RN  RRRRY  CONTRININO  THE 

C INPUT  ORTR.  RNO  THE  OUTPUT  13  TO  8E  OETERfllNEO . IF  I3I0N  13  -1  THE  X 


C VECTOR  IS  THE  OUTPUT  FROfl  THE  NETWORK.  RNO  THE  INPUT  IS  TO  BE  OETERfllNEO. 

C N 13  THE  NUHBER  OF  ORTR  POINTS  RNO  MUST  BE  EOURL  TO  2»»K  WHERE  K 13  RN 

C INTEOER  ORERTER  THAN  0 (RPPENO  ZEROS  TO  THE  ORTR  IF  NECE3SRRY . ) • X MUST 

C BE  OIHEN3IONEO  RT  LEAST  2»N  IN  THE  HRIN  PROORRH.  THIS  ROUTINE  USES  THE 

C COOLEY-TUKEY  FR3T  FOURIER  TRANSFORM  (FOURi)f  HOWEVER.  IT  13  NOT  NECESSARY 
C TO  PUT  THE  ORTR  INTO  COMPLEX  FORM  RS  REOUIREO  8Y  FOUR1 . THE  INPUT  ORTR  IB 

C R33UME0  TO  BE  IN  THE  LEFT  HALF  OF  THE  X ARRAY.  RNO  THE  TRRN3F0RME0  ORTR  18 

C RETURNEO  TO  THE  LEFT  HALF  OF  THE  X RRRRY.  REPLRCINO  THE  INPUT.  IT  13  NOT 
C NECESSARY  TO  00  RNYTHINO  WITH  THE  RIOHT  HALF  OF  THE  X RRRRY.  OT  IS  THE 

C 3RMPILIN0  RATE  OF  THE  ORTR  (SEC).  BF  IS  THE  BRERK  FREQUENCY  (HZ).  CF 

C 13  RN  OPTIONAL  SHARP  CUT  OFF  FREQUENCY  (HZ). 

C 

0IMEN3IQN  XU) 

N0=2»N 

00  1 1 = 1 .N 

X(N0-2»IM)=X(N-I*l) 

1 X(N0-2«I*2)=0.0 
CALL  F0UR1 ( X • M • 1 ) 

FN=FLORT(N) 

I2=N/2*1 

1=2 

2 RMPR=0 .0 
RHP  1=0  »0 

W=FLOAT(I-1)/(FN»OT) 

IF(W.OT.CF)  00  TO  6 
IF(ISION)  4.4.3 

3 RHPR=1 .0/(1 .0»( W/BF )»»2 ) 

RMPI=-(W/BF)/( 1.0*(W/8F)»»2) 

00  TO  6 

4 AMPR=l.O 
AHPI=W/BF 

6 TEMPI=X( 2»I ) 

X(2*I)=X(2»l)«RHPR*X(2»I-l)»RHPI 
X(  2» I- 1 )=X( 2» 1-1 )»RMPR-TEMPI "RMPI 
IFd.E0.I2)  00  TO  6 
X(N0-2»I*4)=-X(2»I) 

X(N0-2»I*3)=X(2»I-l) 

! = !♦! 

00  TO  2 

6 CALL  FOURl(X.N.-l) 

00  7 1=1. N 

7 X( I )=X( 2«I-1 )/FN 
RETURN 

ENO 
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SUBROUTINE  RCN2CX.N .OT.BF .CF.ISION) 

C 

C A SIMPLE  HI0H-PA3S  RC  NETWORK  RS  8H0WN  BELOMi 


C 

C 

C C 

C 0 \\ 0 FREQUENCY  TRANSFER 

C \ FUNCTION 

C \ 

C INPUT  X OUTPUT  FC JW)=(  JW/(8F*JW) )»»I5ION 

C X R 

C X WHERE i 

C \ BF=1/(R»C) 

C 0 0 J=3QRT( -1 ) 

C WrFREQUENCY  (HZ) 

C 

C 

C I 3 ION  IS  *1  OR  -It  IF  ISION  IS  ♦ ! THE  X VECTOR  IS  AN  ARRAY  CONTAININO  THE 


C INPUT  ORTA.  RNO  THE  OUTPUT  IS  TO  SE  DETERMINED.  IF  ISION  IS  -l  THE  X 

C VECTOR  IS  THE  OUTPUT  FROM  THE  NETWORK.  ANO  THE  INPUT  IS  TO  BE  DETERMINED. 

C N IS  THE  NUMBER  OF  ORTA  POINTS  RNO  MUST  BE  EQUAL  TO  2»»K  WHERE  K 19  AN 

C INTEOER  ORERTER  THAN  0 (RPPENO  ZEROS  TO  THE  ORTR  IF  NECESSARY.).  X MUST 

C BE  OIMENSIONEO  AT  LEAST  2-N  IN  THE  MAIN  PROORAM.  THIS  ROUTINE  USES  THE 

C COOLEY-TUKEY  FAST  FOURIER  TRANSFORM  (FOURl)i  HOWEVER.  IT  IS  NOT  NECESSARY 
C TO  PUT  THE  ORTA  INTO  COMPLEX  FORM  AS  REQUIREO  BY  FOUR1.  THE  INPUT  ORTA  18 

C ASSUMED  TO  BE  IN  THE  LEFT  HALF  OF  THE  X ARRAY.  ANO  THE  TRANSFORMED  ORTA  IS 

C RETURNED  TO  THE  LEFT  HALF  OF  THE  X ARRAY,  REPLACINO  THE  INPUT.  IT  IS  NOT 
C NECESSARY  TO  00  ANYTHINO  WITH  THE  RIOHT  HALF  OF  THE  X ARRAY.  OT  IS  THE 

C SAMPILINO  RATE  OF  THE  ORTA  (SEC).  BF  IS  THE  BREAK  FREQUENCY  (HZ).  CF 

C IS  AN  OPTIONAL  SHARP  CUT  OFF  FREQUENCY  (HZ). 

C 

DIMENSION  XU) 

N0=2»N 

00  1 1=1  .N 

X(N0-2»I*1)=X(N-I*i) 

1 X(N0-2» I ♦2)=0 .0 
CALL  FOURKX.N.l) 

FN=FLOAT(N) 

I2=N/2*1 

1=2 

2 AMPR=0.0 
AMPI=0.0 

WrFLOAT(I-l)/(FN»OT) 

IF(W.OT.CF)  00  TO  6 
IF(ISION)  4.3.3 

3 AHPR=l .0/( 1 *0+( 6F/U )a»2) 

AMPI=1.0/(BF/W*W/BF) 

00  TO  6 

4 AHPRrl.O 
AMPI=-8F/W 

6 TEMPI=X( 2» I ) 

X(2»I)=X(2»I)»AHPR^X(2»I-1)»AMPI 
X(2«I-l)=X(2»I-l)»AMPR-TEMPI»AMPI 
IFd.EQ.I2)  00  TO  8 
X( NO-2* 1*4 )s-X( 2=1 ) 

X(N0-2»I*3)=X(2»l-l) 

1=1*1 
00  TO  2 

6 CALL  FOURl(X.N.-l) 

00  7 1=1. N 

7 X( I )*X( 2»I-i )/FN 
RETURN 

ENO 
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SUBROUTINE  RCN3< X .N .OT .BF .A .CF . ISION ) 


r 


\ 


c 

C RN  RC  NETWORK  RS  SHOWN  BELOWi 
C 
C 
C 

C 0 

C 
C 
C 
C 
C 

C INPUT 

C 

c 
c 

C 0 

c 

c 

C IS  I ON  IS  M OR  -li  IF  I SION  IS  *\  THE  X VECTOR  IS  RN  RRRRT  CONTAINING  THE 
C INPUT  DATA.  RND  THE  OUTPUT  IS  TO  8E  DETERMINED . IF  I3I0N  IS  -1  THE  X 

C VECTOR  IS  THE  OUTPUT  FROM  THE  NETWORK.  RNO  THE  INPUT  IS  TO  BE  OETERMINEO. 

C N IS  THE  NUMBER  OF  ORTR  POINTS  RNO  MUST  BE  EQUAL  TO  2»»K  WHERE  K IS  RN 

C INTEGER  GREATER  THAN  0 (APPEND  ZEROS  TO  THE  ORTR  IF  NECESSARY . ) . X MUST 

C BE  OIMENSIONEO  RT  LEAST  2»N  IN  THE  MAIN  PROORRM.  THIS  ROUTINE  USES  THE 

C COOLET-TUKEY  FAST  FOURIER  TRANSFORM  ( FOUR  1 ) f HOWEVER.  IT  IS  NOT  NECESSARY 
C TO  PUT  THE  ORTA  INTO  COMPLEX  FORM  AS  REQUIRED  BY  FOUR!.  THE  INPUT  ORTA  IS 

C ASSUMEO  TO  BE  IN  THE  LEFT  HALF  OF  THE  X ARRAY,  RNO  THE  TRANSFORMED  ORTA  IS 

C RETURNEO  TO  THE  LEFT  HALF  OF  THE  X ARRAY.  REPLACING  THE  INPUT.  IT  IS  NOT 

C NECESSARY  TO  00  ANYTHINO  WITH  THE  RIOHT  HALF  OF  THE  X ARRRY . OT  IS  THE 

C SAMP  I L I NO  RATE  OF  THE  DATA  (SEC).  CF  IS  AN  OPTIONRL  SHARP  CUT  OFF 

C FREQUENCY  (HZ). 

C 

DIMENSION  XU) 

N0=2»N 

00  i Is l . N 

X(N0-2«IM)sX(N-I*l) 

1 X( NO-2* I+2)=0.0 
CALL  FOURl(X.N.l) 

FN=FLOAT(N) 

I2=N/2M 

1=2 

2 AMPR=0.0 
AMP  1:0*0 

W=FLOAT( 1-1 )/(FN*OT) 

IF(W.OT.CF)  GO  TO  S 
IF( IS  ION)  4.3.3 

3 AHPR=( X .0»A»(W/BF)«»2)/(  1 .0*(  A-W/BF )»«2 ) 

AMP  I sW»BF» ( 1 .0-A)/(BF»»2»(A»W)»»2) 

00  TO  6 

4 AMPR=( 1 . 0*A»( W/BF )»»2 )/( 1 .0*( W/BF  )«»2 ) 

AMPI=(A-l.O)/(BF/W*W/BF) 

6 TEMPI=X( 2» I ) 

X(2»I  )=Xf  2" I )»AMPR»X( 2» I- 1 )»AMPI 
X( 2* I - l ) = X( 2» I - 1 )»AMPR-TEMPI»AMPI 
IFd.EQ.I2)  00  TO  6 
X(N0-2«U4)s-X(2»I) 

X(N0-2»I*3)=X<2»I-1) 

1 = I ♦ l 
00  TO  2 

6 CALL  FOURl(X.N.-l) 

00  7 Isl.N 

7 X( I )=X( 2* I-i )/FN 
RETURN 

END 


R1 

VVVV — 0 

\ 

\ 

X F(JW)s(( 

X R2 
X 

\ OUTPUT 

---  C 

\ 

♦ o 


FREQUENCY  TRANSFER 
FUNCTION 

1/A)«(BF*JW)/(BF/A*JW))«»!3I0N 

WHERE i 
A=ioRt/R2 
BFsl/( R2»C ) 

J=SQRT(-1) 

W=FREQUENCY  (HZ) 
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